

function preproc3_ica(sb)


gavgdir = '/home/predatt/anatod/ISI/MEG/GAVG';
cd(gavgdir)
load subjname
load dataset

rootdir = '/home/predatt/anatod/ISI';
datadir = '/home/predatt/anatod/ISI/raw'; 

%% do ICA

infofile = strcat(subjname{sb},'_datainfo.mat'); 
subdir = fullfile(rootdir,'MEG',subjname{sb}); 
cd(subdir);
load(infofile)
cd(datadir)


    % preprocess data
    
       
        cfg = [];
        cfg = datainfo.basicinfo;
        cfg.trl = datainfo.trl;
        
        cfg.detrend ='no';
        cfg.continuous = 'yes';
        cfg.channel = {'MEG', 'MEGREF'};
    
        data = ft_preprocessing(cfg);

    % denoise using third order gradiometer
        cfgdenoise.gradient = 'G3BR';
        data = ft_denoise_synthetic(cfgdenoise, data);

        if sb == 17; % subject had a bad sensor, interpolate before going on

            cfg = [];
            cfg = datainfo.basicinfo;
            cfg.trl = datainfo.trl;
            cfg.method = 'template';
            cfg.layout = 'CTF275.lay';
            cfg.neighbours = ft_prepare_neighbours(cfg);
            cfg.missingchannel = {'MRF32'};
            cfg.method = 'nearest';
            [data] = ft_channelrepair(cfg, data);

        end

    % resample to speed up the ICA decomposition
    cfg            = [];
    cfg.detrend    = 'no';
    cfg.resamplefs = 150;
    data           = ft_resampledata(cfg, data);
   
    % perform the independent component analysis, i.e. decompose the downsampled original data
    % this results in a mixing matrix and in the component timecourses
    % (i.e. the umnixing matrix)
    cfg            = [];
    cfg.channel    = {'MEG'};
    cfg.numcomponent = 60;
    comp           = ft_componentanalysis(cfg, data);

    datainfo.comp.topo = comp.unmixing;
    datainfo.comp.topolabel = comp.topolabel;
    datainfo.comp.do_ica = 'yes';
    
    cd(subdir)
    icafile = strcat(subjname{sb},'_ICA.mat'); 
    save(icafile,'comp');

    clear comp;
     cd(datadir)
    
    %%%% ECG %%%%
    cfg                      = [];
    cfg = datainfo.basicinfo;
    cfg.trl = datainfo.trl;
    cfg.continuous           = 'yes';
    cfg.artfctdef.ecg.pretim = 0.25; %
    cfg.artfctdef.ecg.psttim = 0.5; %
    cfg.artfctdef.ecg.method  = 'zvalue'; % peak-detection method
    cfg.artfctdef.ecg.cutoff  = 3; % peak-threshold
    cfg.artfctdef.ecg.channel = 'EEG058'; % channel EEG058 was the ECG channel in this study
    cfg.artfctdef.ecg.feedback = 'no';
    [cfg, artifact]          = ft_artifact_ecg(cfg); %without feedback

    % preproces the data around the QRS-complex, i.e. read the segments of raw data containing the ECG artifact
    cfg             = [];
    cfg.dataset     = datainfo.basicinfo.dataset;
    cfg.continuous  = 'yes';
    cfg.padding     = 10;
    cfg.dftfilter   = 'yes';
    cfg.demean      = 'yes';
    cfg.trl         = [artifact zeros(size(artifact,1),1)]; %t=0 corresponds to the beginning of the Q (small negative dip before the big positive R)
    if size(cfg.trl,1)>500 cfg.trl = cfg.trl(1:500,:); end; % 500 heartbeats is plenty for identifying the component
    cfg.channel     = {'MEG' 'EEG058'};%
    
    data_ecg      = ft_preprocessing(cfg);

    % resample to speed up the decomposition and frequency analysis, especially useful for 1200Hz MEG data
    cfg            = [];
    cfg.detrend    = 'no';
    cfg.resamplefs = 150;
    data_ecg       = ft_resampledata(cfg, data_ecg);

    % decompose the ECG-locked datasegments into components, using the previously found unmixing matrix
    cfg           = [];
    cfg.unmixing  = datainfo.comp.topo;
    cfg.topolabel = datainfo.comp.topolabel;
    comp_ecg      = ft_componentanalysis(cfg, data_ecg);

    % average the components timelocked to the QRS-complex
    cfg      = [];
    timelock_ecg = ft_timelockanalysis(cfg, comp_ecg);

    for cmp=1:length(comp_ecg.label) % go through all Nchan components and calculate their maximumum amplitude
        ampl_ecg(cmp) = max(abs(timelock_ecg.avg(cmp,:)));
    end

    [ampl_ecg ampl_ecg_i] = sort(ampl_ecg,'descend'); %sort them from high to low
    datainfo.ecg.timelock = timelock_ecg; datainfo.ecg.ampl = ampl_ecg; datainfo.ecg.ampl_i = ampl_ecg_i;
 
   cd(datadir)
   
    %%%% EOG %%%%
    % go back to the raw data on disk and detect the peaks in the VEOG channel
    cfg                       = [];
    cfg.trl                   = datainfo.trl;
    cfg.dataset               = datainfo.basicinfo.dataset;
    cfg.continuous  = 'yes';
    cfg.interactive           = 'yes';
    cfg.artfctdef.eog.channel = {'EEG057'};
    cfg.artfctdef.eog.cutoff  = 4; %peak-threshold
    [cfg, artifact]           = ft_artifact_eog(cfg);

    % preproces the data around the eye-blink
    cfg           = [];
    cfg.dataset   = datainfo.basicinfo.dataset;
    cfg.continuous  = 'yes';
    cfg.padding   = 10;
    cfg.dftfilter = 'yes';
    cfg.demean       = 'yes';
    cfg.trl       = [artifact zeros(size(artifact,1),1)];
    if size(cfg.trl,1)>500 cfg.trl = cfg.trl(1:500,:); end; % 500 eye blinks is plenty for identifying the component
    cfg.channel   = {'MEG' 'EEG057'};
    data_eog      = ft_preprocessing(cfg);

    % resample to speed up the decomposition and frequency analysis, especially useful for 1200Hz MEG data
    cfg            = [];
    cfg.detrend    = 'no';
    cfg.resamplefs = 150;
    data_eog       = ft_resampledata(cfg, data_eog);

    % decompose the VEOG-locked datasegments into components, using the previously found unmixing matrix
    cfg           = [];
    cfg.unmixing  = datainfo.comp.topo;
    cfg.topolabel = datainfo.comp.topolabel;
    comp_eog      = ft_componentanalysis(cfg, data_eog);
    
    % average the components timelocked to the eye-blink
    cfg      = [];
    cfg.vartrllength =2;
    timelock_eog = ft_timelockanalysis(cfg, comp_eog);

    for cmp=1:length(comp_eog.label) % go through all components
        ampl_eog(cmp) = max(abs(timelock_eog.avg(cmp,:)));
    end

    
    [ampl_eog ampl_eog_i] = sort(ampl_eog,'descend');
    datainfo.eog.timelock = timelock_eog; datainfo.eog.ampl = ampl_eog; datainfo.eog.ampl_i = ampl_eog_i;
   
    
 cd(subdir)
 save(infofile,'datainfo');

 
end
